# example of raw sattelite processing 
options(warn=-1)
library("raster")
library("sf")
library("tidyverse")
library("rasterVis")

config.ndwi_treshold = 0 # allowed: -0.05, 0, 0.1 by article 
config.mask <- st_read("data/mask/tile.geojson")

b3 <- raster("data/sattelite/2023/03-20/2023-03-20-00:00_2023-03-20-23:59_Sentinel-2_L2A_B03_(Raw).tiff")
b4 <- raster("data/sattelite/2023/03-20/2023-03-20-00:00_2023-03-20-23:59_Sentinel-2_L2A_B04_(Raw).tiff")
b5 <- raster("data/sattelite/2023/03-20/2023-03-20-00:00_2023-03-20-23:59_Sentinel-2_L2A_B05_(Raw).tiff")
b8 <- raster("data/sattelite/2023/03-20/2023-03-20-00:00_2023-03-20-23:59_Sentinel-2_L2A_B08_(Raw).tiff")

b3 <- raster::mask(b3, config.mask)
b4 <- raster::mask(b4, config.mask)
b5 <- raster::mask(b5, config.mask)
b8 <- raster::mask(b8, config.mask)

ndwi <- (b3 - b8)/(b3 + b8)
# make binary reclassification to 2 classes: from -1 to 0.3 => 0, from 0.3 to 1 => 1
cndwi <- reclassify(ndwi, c(-1, config.ndwi_treshold, 0, config.ndwi_treshold, 1, 1))

mask_ndwi <- rasterToPolygons(cndwi, fun = function(x){x == 1}, digits=6, dissolve = T)

b3 <- raster::mask(b3, mask_ndwi)
b4 <- raster::mask(b4, mask_ndwi)
b5 <- raster::mask(b5, mask_ndwi)
b8 <- raster::mask(b8, mask_ndwi)


# calc chlorophyll index (normalized)
ndci <- (b5-b4)/(b5 + b4)

# plot output NDCI
pal <- terrain.colors(4)
plot(ndci, breaks = c(0, 0.2, 0.4, 1), col=pal)
# process values
vals <- as.vector(na.omit(ndci@data@values))
stats.ndci = list(
  mean = mean(vals),
  sd = sd(vals),
  stderr = sd(vals)/sqrt(length(vals)),
  min = min(vals),
  max = max(vals)
)

hist(vals)

# NDCI to ug/l, exact as mg/m^3
# ndci to microg/l by Xing Wang et al., 2018
# y = 20.877 * ncdi + 3.273
ugl <- 20.877*ndci + 3.273
pal <- terrain.colors(5)
ugl@data@values[ugl@data@values < 0] <- 0
plot(ugl, col=terrain.colors(15))
vals <- as.vector(na.omit(ugl@data@values))
vals[vals < 0] <- 0
stats.chl = list(
  mean = mean(vals),
  sd = sd(vals),
  stderr = sd(vals)/sqrt(length(vals)),
  min = min(vals),
  max = max(vals)
)

levelplot(ugl, main="Хлорофилл-а, 03-10", contour = F, margin = F, col.regions = terrain.colors(100))

